用BayesSpace提高分群和基因分辨率
作者 | 周运来
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树。
生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。
最近开发的空间基因表达技术，如空间转录组学可以在保留空间背景的同时，全面测量转录组谱。然而，现有的空间基因表达数据分析方法往往不能有效利用空间信息，不能解决技术分辨率有限的问题。
在空间数据的分析中，痛点像单细胞转录组一样也是：聚类。我们希望聚出来的类不仅在uamp、tsne这样的空间中是清楚明了的，也希望在空间坐标中也是这样。同时也希望找到空间特异的基因表达模式。今天给大家介绍的BayesSpace，正是在做这方面的努力。
BayesSpace是一种完全贝叶斯统计方法，用于聚类分析来增强空间转录组数据的分辨率，并且可以无缝集成到当前的转录组分析工作流程中。BayesSpace利用多元t分布对基因表达的低维表示进行建模，然后通过Potts模型合并空间信息，这使得相邻点属于同一簇。这种方法借鉴了用于图像分析的成熟的空间统计方法。有效地利用空间信息可以极大地提高了空间数据聚类的性能。
在算法考虑并鼓励邻居点属于同一个簇的原则也可以用来提高点的分辨率到“子点（sub-spot）”水平。这里鼓励子点属于同一个簇。子点水平表达有时候并不能直接观察到，但可以通过反复抖动点水平表达值，通过MCMC进行估计。在实践中，bayesspace分辨率增强的聚类可以识别重要的生物结构，否则可能会错过。
我们从10X下载他们做的示例数据，并安装好了BayesSpace，下面让我们探索一番，并与Seurat的聚类结果做一个比较。
library(SpatialExperiment)
library(SingleCellExperiment)
library(ggplot2)
library(BayesSpace)
library(Seurat)
library(cowplot)
allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",
"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
"#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
"#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887") # 我最爱的颜色搭配
分别用BayesSpaceh和Seurat读入数据，这里小哥给Seurat的函数都加Seurat::以表示来自Seurat，以此提示亲们函数的来源。同时创建sce和seo两个对象。
sce <- readVisium("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\")
sce
seo<-Seurat::Load10X_Spatial("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\")
seo
先走我们熟悉的Seurat空间流程：
plot1 <- Seurat::VlnPlot(seo, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- Seurat::SpatialFeaturePlot(seo, features = "nCount_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)
Seurat 默认聚类方法：
seo <- Seurat::SCTransform(seo, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
seo <- Seurat::RunPCA(seo, assay = "SCT", verbose = FALSE)
seo <- Seurat::FindNeighbors(seo, reduction = "pca", dims = 1:30)
seo <- Seurat::FindClusters(seo, verbose = FALSE)
seo <- Seurat::RunUMAP(seo, reduction = "pca", dims = 1:30)
p1 <- Seurat::DimPlot(seo, reduction = "umap", label = TRUE)
p2 <- Seurat::SpatialDimPlot(seo, label = TRUE, label.size = 3)
plot_grid(p1, p2)
下面我们来用BayesSpace的一般聚类方法，对数据进行聚类。
set.seed(1314)
?spatialPreprocess
sce <- scater::logNormCounts(sce)
sce <- spatialPreprocess(sce, platform="Visium",
n.PCs=7, n.HVGs=2000, log.normalize=FALSE)
sce <- qTune(sce, qs=seq(2, 10), platform="Visium", d=7)
qPlot(sce)
sce
class: SingleCellExperiment
dim: 33538 3813
metadata(2): BayesSpace.data chain.h5
assays(2): counts logcounts
rownames(33538): MIR1302-2HG FAM138A ... AC213203.1 FAM231C
rowData names(3): gene_id gene_name is.HVG
colnames(3813): CAGGATCCGCCCGACC-1 CACGATTGGTCGTTAA-1 ... CGCAGTTCTATCTTTC-1 GACAGGTAATCCGTGT-1
colData names(9): spot in_tissue ... cluster.init spatial.cluster
reducedDimNames(1): PCA
altExpNames(0):
推荐聚类个数：
虽然看到聚成8个是可以的，但是为了和Seurat有一定程度的比较，我们决定用10。
sce <- spatialCluster(sce, q=10, platform="Visium", d=7,
init.method="mclust", model="t", gamma=2,
nrep=1000, burn.in=100, # nrep建议大于10000 ，奈何我们为了演示只能牺牲点了
save.chain=TRUE)
比较一下二者的区别，为了颜色更加明显，我们用黑色的背景：
p3 <-SpatialDimPlot(seo, label = T, label.size = 3)+ ggtitle("Seurat") +DarkTheme()
p4 <- clusterPlot(sce,label = "spatial.cluster" )+ ggtitle("BayesSpace")+DarkTheme()
plot_grid(p3, p4)
可以看到，Seurat默认的聚类方法在空间上有一定的混杂，虽然BayesSpace的一般聚类也有一些。但是BayesSpace还为我们准备了spatialEnhance()函数，可以增强主成分的分辨率，并将这些pc和预测的子点分辨率的聚类标签添加到新的聚类单元中。正如上面我们对spatialCluster()的演示一样，我们在这个示例中使用的迭代次数(nrep=1000)要少于我们在实践中推荐的次数(nrep=100000或更大)。
sce.enhanced <- spatialEnhance(sce, q=10, platform="Visium", d=7,
model="t", gamma=2,
jitter_prior=0.3, jitter_scale=3.5,
nrep=1000, burn.in=100, # 为了向内存妥协
save.chain=T)
增强的singlecelexperimental包括原始sce中的父点索引(spot.idx)以及子点索引。它向原始spot坐标添加偏移量，并提供增强的集群标签(spatial.cluster)。
DataFrame with 6 rows and 9 columns
spot.idx subspot.idx spot.row spot.col row col imagerow imagecol spatial.cluster
<numeric> <integer> <integer> <integer> <numeric> <numeric> <numeric> <numeric> <numeric>
subspot_1.1 1 1 1 49 1.333333 49.3333 4362.18 10258.5 4
subspot_2.1 2 1 0 50 0.333333 50.3333 4125.18 10395.5 6
subspot_3.1 3 1 1 51 1.333333 51.3333 4363.18 10531.5 4
subspot_4.1 4 1 0 52 0.333333 52.3333 4126.18 10668.5 6
subspot_5.1 5 1 1 53 1.333333 53.3333 4363.18 10804.5 4
subspot_6.1 6 1 0 54 0.333333 54.3333 4126.18 10941.5 4
我们也将两者做一个比较：
clusterPlot(sce.enhanced,label = "spatial.cluster" )+ ggtitle("enhanced")+DarkTheme()
bayesspace可视化基因可以提高基因表达的分辨率。BayesSpace处理基因表达矩阵的主成分，spatialEnhance()函数计算增强分辨率的PC向量。增强的基因表达不是直接计算，而是使用回归算法输入，对于每个基因训练一个使用每个位点的PC载体的模型来预测位点水平的基因表达，并使用拟合的模型从子位点PCs中预测子点表达。
基因表达增强是在enhanceFeatures()函数中实现的。BayesSpace默认情况下使用xgboost预测表达式，但是线性回归和狄利克雷回归也可以通过model参数使用。在使用xgboost时，我们建议通过将nrounds参数设置为0来自动调优它，尽管这是以增加运行时为代价的(实际中比预先指定的nrounds要慢约4倍)。
enhanceFeatures()可用于计算所有基因或感兴趣基因子集的子点水平表达。在这里，我们将通过增强四种标记基因的表达来证明:PMEL(黑色素瘤)，CD2 (t细胞)，CD19 (b细胞)和COL1A1(成纤维细胞)。
markers <- c("PMEL", "CD2", "CD19", "COL1A1")
sce.enhanced <- enhanceFeatures(sce.enhanced, sce,
feature_names=markers,
nrounds=0)
enhanced.plots <- purrr::map(markers, function(x) featurePlot(sce.enhanced, x))
spot.plots <- purrr::map(markers, function(x) featurePlot(sce, x))
patchwork::wrap_plots(c(enhanced.plots, spot.plots), ncol=4)
虽然我们的marker只是根据推测来选的，很明显可以看出，enhanced之后基因的表达模式更加凸显，而未处理的基因表达规律不强。
References
[1] BayesSpace enables the robust characterization of spatial gene expression architecture in tissue sections at increased resolution: https://www.biorxiv.org/content/10.1101/2020.09.04.283812v1[2] https://research.fredhutch.org/gottardo/en/blog/bayesspace.html
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注